Week12 instructor

   

Week 12: Multi-band Raster

In [10]:

import os
import rasterio
import numpy as np
import rasterio.mask
import rasterio.plot
import geopandas as gpd
import matplotlib.pyplot as plt

import os import rasterio import numpy as np import rasterio.mask import rasterio.plot import geopandas as gpd import matplotlib.pyplot as plt

Reading Multiple Band¶

In [11]:

os.getcwd()

os.getcwd()

Out[11]:

'e:\\Clark\\Spring2025\\Github_Courses\\geog386_pythonProgramming\\docs\\Lectures\\Week12_multiple_bands'

In [12]:

dem_path = os.path.join(os.getcwd(),  "landsat.tif")
dem_path

dem_path = os.path.join(os.getcwd(), "landsat.tif") dem_path

Out[12]:

'e:\\Clark\\Spring2025\\Github_Courses\\geog386_pythonProgramming\\docs\\Lectures\\Week12_multiple_bands\\landsat.tif'

In [13]:

landsat = rasterio.open(dem_path)

landsat = rasterio.open(dem_path)

Check metadata of raster¶

In [14]:

landsat.meta

landsat.meta

Out[14]:

{'driver': 'GTiff',
 'dtype': 'int16',
 'nodata': -32768.0,
 'width': 776,
 'height': 776,
 'count': 7,
 'crs': CRS.from_wkt('PROJCS["WGS 84 / UTM zone 17N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32617"]]'),
 'transform': Affine(30.0, 0.0, 202215.0,
        0.0, -30.0, 4031715.0)}

Visualizing Multiple Bands¶

Name	Wavelength (µm)	Description
Band1: Coastal	0.44 - 0.45	Coastal and aerosol
Band 2: Blue	0.45 - 0.51	Water body analysis
Band 3: Green	0.53 - 0.59	Vegetation and soil discrimination
Band 4: Red	0.64 - 0.67	Vegetation analysis
Band 5: Near Infrared (NIR)	0.86 - 0.88	Biomass content and water body analysis
Band 6: Shortwave Infrared (SWIR) 1	1.6 - 1.7	Soil moisture, vegetation, and snow differentiation
Band 7: Shortwave Infrared (SWIR) 2	2.2 - 2.3	Penetrates clouds, smoke, and haze

Read a specific band¶

What is Numpy array¶

• Image Source

Raster image includes many different bands. Each band can be represent as array in Python

• Data Source

In [15]:

rasterio.plot.show((landsat,5),cmap="Greys_r")

rasterio.plot.show((landsat,5),cmap="Greys_r")

Out[15]:

<Axes: >

rasterio.plot.show(landsat)

Read all bands as 3D¶

In [16]:

# Read all bands
img_array = landsat.read()
img_array.shape

# Read all bands img_array = landsat.read() img_array.shape

Out[16]:

(7, 776, 776)

Read single band as 2D¶

In [17]:

img_array = landsat.read(2)
img_array.shape

img_array = landsat.read(2) img_array.shape

Out[17]:

(776, 776)

Calculate the band statistics¶

In [18]:

# Read all bands
bands = landsat.read()

# Calculate statistics for each band
stats = []
for band in bands:
    stats.append({
        'min': band.min(),
        'mean': band.mean(),
        'median': np.median(band),
        'max': band.max()})

# Show stats for each channel
stats

# Read all bands bands = landsat.read() # Calculate statistics for each band stats = [] for band in bands: stats.append({ 'min': band.min(), 'mean': band.mean(), 'median': np.median(band), 'max': band.max()}) # Show stats for each channel stats

Out[18]:

[{'min': np.int16(-32768),
  'mean': np.float64(6966.947397438623),
  'median': np.float64(7781.0),
  'max': np.int16(25336)},
 {'min': np.int16(-32768),
  'mean': np.float64(7084.704948387182),
  'median': np.float64(7862.0),
  'max': np.int16(27531)},
 {'min': np.int16(-32768),
  'mean': np.float64(7811.792555664789),
  'median': np.float64(8603.0),
  'max': np.int16(30420)},
 {'min': np.int16(-32768),
  'mean': np.float64(7482.2932813662455),
  'median': np.float64(8191.0),
  'max': np.int16(32261)},
 {'min': np.int16(0),
  'mean': np.float64(16486.904418309598),
  'median': np.float64(16946.0),
  'max': np.int16(32767)},
 {'min': np.int16(0),
  'mean': np.float64(12657.95134312892),
  'median': np.float64(12898.0),
  'max': np.int16(32767)},
 {'min': np.int16(-32768),
  'mean': np.float64(9008.560709825699),
  'median': np.float64(9696.0),
  'max': np.int16(32532)}]

Visualizing single band¶

In [29]:

band_names = ["Coastal", "Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"]

band_names = ["Coastal", "Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"]

In [20]:

fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(8, 10))
axes = axes.flatten()  # Flatten the 2D array of axes to 1D for easy iteration

for band in range(1, landsat.count):
    data = landsat.read(band)
    ax = axes[band - 1]
    im = ax.imshow(data,cmap="gray")
    ax.set_title(f"Band {band_names[band - 1]}")
    fig.colorbar(im, ax=ax, label="Reflectance", shrink=0.5)

plt.tight_layout()
plt.show()

fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(8, 10)) axes = axes.flatten() # Flatten the 2D array of axes to 1D for easy iteration for band in range(1, landsat.count): data = landsat.read(band) ax = axes[band - 1] im = ax.imshow(data,cmap="gray") ax.set_title(f"Band {band_names[band - 1]}") fig.colorbar(im, ax=ax, label="Reflectance", shrink=0.5) plt.tight_layout() plt.show()

Stacking raster bands¶

Band Combinations for Landsat 8

• The bands are arranged in the Red-Green-Blue order

In [21]:

red_band = landsat.read(6)
green_band = landsat.read(5)
blue_band = landsat.read(2)

red_band = landsat.read(6) green_band = landsat.read(5) blue_band = landsat.read(2)

In [22]:

def normalize(array):
    array_max = array.max()
    return (array/array_max)

n_blue = normalize(blue_band)
n_red = normalize(red_band)
n_green = normalize(green_band)

# Stack the bands into a single array
rgb = np.dstack((n_red, n_green, n_blue)).clip(0, 1)

# Plot the stacked array
plt.figure(figsize=(8, 8))
plt.imshow(rgb)
plt.title("Band Combination: RGB")
plt.show()

def normalize(array): array_max = array.max() return (array/array_max) n_blue = normalize(blue_band) n_red = normalize(red_band) n_green = normalize(green_band) # Stack the bands into a single array rgb = np.dstack((n_red, n_green, n_blue)).clip(0, 1) # Plot the stacked array plt.figure(figsize=(8, 8)) plt.imshow(rgb) plt.title("Band Combination: RGB") plt.show()

In [23]:

rgb.shape

rgb.shape

Out[23]:

(776, 776, 3)

Band Math (NDVI Calculation)¶

NDVI is calculated as:

NDVI = (NIR - Red) / (NIR + Red)

In [24]:

nir = landsat.read(5)
red_band = landsat.read(4)

ndvi = (nir - red_band) / (nir + red_band)

plt.figure(figsize=(8, 8))
plt.imshow(ndvi, cmap="RdYlGn", vmin=-1, vmax=1)
plt.colorbar(label="NDVI", shrink=0.5)
plt.title("NDVI")
plt.xlabel("Column #")
plt.ylabel("Row #")
plt.show()

nir = landsat.read(5) red_band = landsat.read(4) ndvi = (nir - red_band) / (nir + red_band) plt.figure(figsize=(8, 8)) plt.imshow(ndvi, cmap="RdYlGn", vmin=-1, vmax=1) plt.colorbar(label="NDVI", shrink=0.5) plt.title("NDVI") plt.xlabel("Column #") plt.ylabel("Row #") plt.show()

Extract vegetation area¶

In [25]:

ndvi_clean = np.where(ndvi > 0.2, ndvi, np.nan)
# ndvi_clean.plot(cmap="Greens", vmin=0.2, vmax=0.8)

nd_image = plt.imshow(ndvi_clean, cmap="YlGn", vmin=0.2, vmax=1)

cbar = plt.colorbar(nd_image, shrink = 0.5)
cbar.set_label("NDVI")

plt.title("Vegetated Areas")
plt.xlabel("Column #")
plt.ylabel("Row #")
plt.show()

ndvi_clean = np.where(ndvi > 0.2, ndvi, np.nan) # ndvi_clean.plot(cmap="Greens", vmin=0.2, vmax=0.8) nd_image = plt.imshow(ndvi_clean, cmap="YlGn", vmin=0.2, vmax=1) cbar = plt.colorbar(nd_image, shrink = 0.5) cbar.set_label("NDVI") plt.title("Vegetated Areas") plt.xlabel("Column #") plt.ylabel("Row #") plt.show()

Visualizing Contours on NDVI¶

Overlay contour lines on an NDVI heat map to highlight zones of similar NDVI values. Contours can help identify areas with different vegetation densities and enhance visual interpretation.

Useful for identifying zones with similar vegetation density, making it easier to delineate high versus low vegetation areas

In [26]:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D

In [27]:

plt.figure(figsize=(10, 8))
plt.imshow(ndvi, cmap="RdYlGn", vmin=-1, vmax=1)
plt.colorbar(label="NDVI", shrink = 0.5)
plt.contour(ndvi, levels=10, colors='black', alpha=0.5)
plt.title("NDVI with Contour Overlay")
plt.show()

plt.figure(figsize=(10, 8)) plt.imshow(ndvi, cmap="RdYlGn", vmin=-1, vmax=1) plt.colorbar(label="NDVI", shrink = 0.5) plt.contour(ndvi, levels=10, colors='black', alpha=0.5) plt.title("NDVI with Contour Overlay") plt.show()

Show the NDVI range and frequency through histogram¶

In [28]:

from rasterio.plot import show_hist
show_hist(
    ndvi, bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", label = "NDVI")

from rasterio.plot import show_hist show_hist( ndvi, bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled', title="Histogram", label = "NDVI")